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BAYESIAN METHODS FOR FLOW PARAMETER ESTIMATES IN 
MAGNETIC RESONANCE IMAGING 



CROSS REFERENCES TO RELATED APPLICATIONS 

This application claims the benefit of United States Provisional Application Serial No. 
60/181,823 filed on February 1 1, 2000. 

FIELD OF THE INVENTION 

This invention relates to the field of signal processing, and more particularly to magnetic 

resonance imaging systems. 

BACKGROUND OF THE INVENTION 

Magnetic Resonance Imaging (MRI) has proven to be a revolutionary diagnostic 
radiological tool, due to its high spatial resolution and excellent discrimination of soft tissues. 
Magnetic Resonance Imaging provides rich information about anatomical structure, enabling 
quantitative anatomical studies of diseases, the derivation of computerized anatomical atlases, as 
well as three-dimensional visualization of internal anatomy, for use in pre-operative and intra- 
operative visualization, and in the guidance of therapeutic intervention. 

The role of dynamic imaging in medicine is of growing importance as a tool for both 
diagnosis and research. Most importantly is the role of flow studies inside the vascular tree: 
these range from flow studies in the heart, major arteries, and through replacement devices like 
valves, to more delicate studies of flow patterns in the brain. 

For very large superficial vessels like the carotid, methods such as Doppler have proved 
very useful. Flow studies of smaller, deeper vessels have necessitated more complicated and 
invasive techniques such as the scanning of exogenously introduced radio labeled substances. 
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Valuable as these techniques have been, the results have lacked resolution to provide detailed 
flow parameters on a small scale. 

Therefore, there is a need to be able to provide higher resolution and detail flow 
parameters for smaller/deeper vessels. 

SUMMARY OF THE INVENTION 

In accordance with the present invention, there is provided a method and system for flow 
parameter estimates in magnetic resonance imaging comprising the following steps: accessing 
magnetic resonance imaging data; inputting a magnetic resonance imaging model function; and, 
using conditional probabilities based on Bayes' Theorem. 

BRIEF DESCRIPTION OF THE DRAWINGS 

A more complete understanding of the present invention may be obtained from 
consideration of the following description in conjunction with the drawings in which: 
Fig. 1 is an overview of a Magnetic Resonance Imaging system; and, 
Fig. 2 is a high level flow chart of the system methodology. 

DETAILED DESCRIPTION OF VARIOUS ILLUSTRATIVE EMBODIMENTS 

A conventional Magnetic Resonance Imaging system is generally depicted in Fig. 1. The 
MRI machine is placed in an RF shielded room 100 where data on the patient will be taken. This 
is to avoid spurious signals from the surroundings: both from external electronics as well as from 
signals generated by the driving electronics of the system itself. The probe part of the MRI 
machine houses magnetic cryostats 102 for keeping the magnetic coils cool and in a peak 
operating range. There are essentially three sets of coils: gradient coils 104, receive coils 106 
and transmit coils 108. The patient is positioned on a patient table 110, which brings the patient 
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into the magnetic field in such a way that the anatomy to be imaged is in the region where the 
field is homogeneous; this region is called the isocenter 1 12. 

A computer system 114, located outside the RF shielded room 100, operates the MRI 
machine. A user requesting specific imaging information via acquisition commands mans the 
computer system 114. The processed data is sent to the viewing section 1 16 on the console 118. 
Here the user can display the images on the monitor 120 and/or can reproduce them on a hard 
copy camera 122. 

To do a targeted scan relevant patient statistics are first entered into the administration 
section on the console 118. The system is then instructed to do the required scan, including the 
geometrical parameters, the imaging method, and the sequence timing. (For new or experimental 
imaging the gradient field strength as a function of time for each direction x, y and z and the RF 
waveform must first be programmed into the memory of the computer system 114. This can also 
be done on the console 118.) The information is then forwarded to an operational area of the 
machine known in general as the "spectrometer." The spectrometer 126 consists of the front- 
end-controller 128 (controlling the magnet, gradients, RF transmitter and receiver, RF coil 
switches, etc.) and the data acquisition around the receiver (the receiver switches and the ADC 
(analog-to-digital converter) 130. 

Before the actual scan begins the spectrometer must perform an initialization sequence: 
tune the synthesizer 132 frequency, o 0 , to the gyromagnetic frequency; tune the RF receiving 
coils; and, the RF power and receiver gain must be adjusted to the specs of the designated 
measurement, etc. 

Now the Spin Echo sequence can start. All components, except for the magnet, turned 
off a selection of the field gradient is. When the field gradient reaches it preset required value 



3 



4683-107 

the RF power amplifier 134 is engaged and its output is switched 136 with frequency co 0 . This 
signal (also known as the excitation pulse signal) is modulated in the waveform generator 138 
and fed into the power amplifier 134. The RF signal, which can be AM or FM modulated, drives 
a current in the transmit coil 108 to produce the required magnetic induction, B RF , within the 
center of the coil. Simultaneously, the receiver coils 106 are detuned and the pre-amplifier 140 
blocked by switch 142 so that the large transmit signal will not burn out the pre-amplifier 140. 
When the RF pulse is of sufficient magnitude, the output line of the power amplifier 134 is 
switched 136 to the matched load again. The output noise of the power amplifier 134 with zero 
input signal is high enough that it could mask the MR signal to be measured and so the switching 
is a necessary function. The selection gradient is now switched off. Because of the high coil 
inductance and the gradient-amplifier voltage the reduction of the gradient field is not 
instantaneous but occurs at some small time interval. As a result of this procedure, spins in a 
thin slice of the target are aligned. 

Next, the dephasing gradient and the phase-encode gradient are switched on and brought 
to the required strength, which has been programmed into the operational specs of the 
spectrometer 126. When the gradient pulses have the required surface value, i.e., the integral of 
the gradient over time, they are switched off and the refocusing pulse can now be applied. This 
180° RF pulse is usually slice selective and can be applied to the same slice as the excitation RF 
pulse. Its waveform is similar to that of the excitation pulse but with twice the amplitude or four 
times the power. In single-slice methods, slice selectivity of the refocusing pulse may be left 
out, making the pulse of shorter duration possible. 

When the refocusing RF pulse cycle is finished the read-out gradient is switched on. 
During this gradient pulse the receiver coil 106 is in a tuned state and the receiver is activated. 
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Concurrently, the magnetization is measured and the ADC 130 samples the received signal in 
predetermined sampling time interval. The system is then allowed a period of time to relax to its 
initial configuration so that the spectrometer 126 can proceed to image the next elemental slice in 
the scan sequence. The process is repeated until the whole scan is completed. Slice sample data 
are sent to the array processor 146 where fast Fourier transforms are performed on the data. 
When the slices are transformed and reassembled, an image results. This image is sent to the 
view section 1 16 on the console 118. 

Usually after this first scan more scans of the same anatomical area follow. These scans 
are taken of slices positioned on the basis of the first image. The user can select the orientation 
and the off-center distance of the slices addressed in such scans. 

The present invention, Bayesian methods for flow parameter estimates magnetic resonance 
imaging is a technique and system, wherein information collected in Nuclear Magnetic 
Resonance (NMR), also known as Magnetic Resonance Imaging (MRI), studies provides provide 
higher resolution and detail. The present invention resolves problems of resolution and deep 
tissue location, thus enabling the use of a non-invasive technique for virtually all flow studies on 

any scale available to MRI. 

The immediate clinical implications are widespread. In the assessment of adequate blood 
flow, such as the determination of cardiac output, flow to the extremities or the head, the present 
invention provides a degree of accuracy and penetration to non-superficial vessels unavailable 
until now. This is especially useful in the evaluation of the extent of thrombotic/embolic strokes 
and myocardial infarctions, both past and in evolution by use of a non-invasive, zero radiation 
dose, real time procedure. 
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In addition to answering the "What is the flow in this vessel?" question, the present 
invention, Bayesian methods for flow parameter, estimates magnetic resonance imaging, enables 
the examination of microscopic flow characteristics on an extremely small scale. This allows the 
examination of flow parameter gradients across a vessel diameter and the quantitative 
determination of turbulence regimes across hydrodynamic impedances such as heart valves, 
plaques, spasms, and malformations such as anneyurisms. This latter capability is very 
important in the design of intravascular prostheses such as valves, vessel replacements, stents, 
etc. 

Nuclear Magnetic Resonance (NMR) is fundamentally a process by which the magnetic 
moments of various nuclei are obtained. The nuclei are first immersed in a fairly strong but 
homogenous magnetic field. This field causes the nuclear magnetic moment of an nucleus to 
precess about an axis parallel to the direction of the field. The frequency with which the 
magnetic moment precesses is called the Larmor frequency. Then, a weak oscillating field is 
applied in an orthogonal direction to the homogenous field. When the frequency of the field 
oscillations is the Larmor frequency, a resonance condition is met and the nuclei are caused to 
jump from one quantum state to another. 

All magnetic resonance experiments involve the transition of the nuclear spins from a 
given state to one of higher energy. In order to be detectable, two conditions must be satisfied. 
First, more nuclei must be in the lower states than in the upper ones, so that the upward 
absorptive transitions outnumber the downward induced emissive transitions which return energy 
to the oscillator. Second, the nuclei in the upper state must have a means of returning to a lower 
state other than by emitting radiation. The first condition is satisfied when the populations of the 
various energy states are in statistical equilibrium at sufficiently low temperatures in a 
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sufficiently strong magnetic field. The second condition is satisfied if there is a sufficiently 
strong coupling between the nuclear moment and the solid or liquid system in which it is 
situated. This coupling leads to nonradiative transitions of the nuclei, the energy being 
transferred to the crystal lattice where it appears as heat. This energy transfer called spin-lattice 
relaxation permits nuclear spins to return toward statistical equilibrium. Another relaxation 
mechanism is called spin-spin relaxation. It is due to the interaction between the magnetic 
moments of two nuclei. The effect of the spin-spin interaction is to cause the transverse 
component — orthogonal to the homogeneous magnetic induction field — of the magnetization to 
relax back to zero. 

The phenomenon of nuclear magnetic resonance has been put to use as a powerful tool 
for discovering new properties of matter. This has been especially true in nuclear physics, 
chemistry, and most recently in medicine where NMR is called MRI. Numerical values of the 
magnetic moments provide information about the structure of the nuclei. 

It is critical to be able to extract and discern all relevant structural information from the 
NMR data in order to have a true and accurate picture of the system of interest. Brethorst 
developed a Bayesian method, based on conditional probabilities, for inverting nuclear magnetic 
resonance measurements-free magnetic induction decay measurements taken on a mixture of 
63% liquid Hydrogen-Deuterium (HD) and Deuterium (D2) at 20.2 K-to obtain critical 
characterization parameters of a target molecular specie. His approach used a Gaussian-like 
Likelihood Function which employed a simple single harmonic frequency Model Function and 
the assumption of background white noise. He then modified the Likelihood Function to include 
multiple harmonic frequencies, assumed not to have any interfering effects within the sample, 
which resulted in a somewhat standard Posterior Probability Function known as the Student-t 
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Distribution. The Bretthorst method accurately estimated quantities such as spin precession 
frequencies, spin decay/relaxation rates, and spin magnetization. His results were markedly 
better than those obtained by stochastic methods. Additionally, the Bretthorst method has the 
ability to resolve very closely spaced frequencies if there is information (evidence) in the data. 
This fine structure is of particular importance in flow studies. 

The use of the Bayesian method has enhanced the information provided by NMR free 
magnetic induction decay measurements on stationary systems. However, it has not been 
modified and applied to dynamic systems or in vivo measurements on patients where 
motion-intentional as well as unintentional-and other than white background noise are critical. 
The present invention, Bayesian methods for flow parameter estimates in magnetic resonance 
imaging, is a technique and system which applies Bayesian techniques to the analysis of complex 
MRI flow data to obtain parameters of interest such as velocity, acceleration, turbulence, and 
phase shifts due to flow gradients. This is done, in part, by uniquely integrating dynamic Model 
Functions into the Bayesian scheme. (Dynamic Model Functions mean parameterized Model 
Functions that are time - dependent in order to specifically address the flow nature of the 
problem.) MRI flow data, in particular in vivo flow data, is infinitely more difficult to resolve 
than in vivo stationary data. The present invention's creative application of Bayesian probability 
methods to in vivo flow data is of extreme benefit in resolving this vastly more complex 
information. 

The Bloch Equations 

If a spin magnet moment m is immersed in a field of magnetic induction B the resulting 
torque on the system is the time derivative of the spin angular momentum L and is given by 
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dL 

dt 



= mxB, 



(2.1) 



Since m = g(e / 2m p c)L = y N L, then 



dm 
dT 



= Y N (mxB). 



(2.2) 



g is the nuclear g-factor whose value depends on the particular nucleus being studied, e is the 
charge of an electron, m p is the rest mass of the spin particle (in the present case, a nucleon — the 
proton / neutron) and c is the speed of light in vacuum. One of the consequences of having the 
magnetic moment proportional to the angular momentum is that when placed in a magnetic field 
it will precess. Indeed, the magnetic moment precesses with a frequency that is dependent upon 
the g-factor or the structure of the nucleus. 

Equation (2.2) can easily be generalized to a system of N nucleons. In this case the 
torque is given by the time derivative of the net magnetic moment. For systems consisting of 
multiple magnetic moments due to individual sources, it is conventional to work with the 
magnetization M, which is the net magnetic moment per unit volume. 



Bloch modified Eq. (2.3) in a phenomenological approach to include spin relaxation of a 
crystalline solid-state system in two ways. First, he included the interactions between all spin 
magnetic moments of the nucleons — the so-called spin-spin interactions. And second, he 
included the relaxation mechanism due to spin-lattice interactions. The resulting Bloch 
equation(s) is 



dM 

dt 



= Yn(MxB). 



(2.3) 



dM 

dt 



= Y N (MxB)- 



M x i + M y j (M z -M 0 ) 



(2.4) 



T 2 T 
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Here M x , M y and M z are the Cartesian components of M, i, j and k are their respective unit 
vectors and M 0 is the equilibrium value of M z , T2 is the spin-spin relaxation time, and Ti 
corresponds to the spin-lattice relaxation time. 

Casting the above in terms of a typical MRI test, the homogenous constant magnetic 
induction field is chosen to be oriented along the positive z-axis, B 0 = B 0 k, of the laboratory 
frame. To this is added an oscillating — rotating — field term having components in the x-y plane, 
Bi = B Ix i + B ly j, i.e. orthogonal to B 0 . This term represents the excitation of the system by an 

RF pulse and will impose its own precession on the magnetization. Finally, a term is added 
which is called the gradient field term. This term has components [J which account for two 
different mechanisms. The excitation pulse that gives rise to Bi is not localized. In MRI 
measurements, excitation of the spins in a selected slice of sample is required. This is 
accomplished by the superposition of a small linear position-dependent field onto the 
homogeneous field. Ideally this would be directed along the z-axis parallel to B 0 . Additional 
components are needed, however, that will account for all the unwanted field deviations. Both of 
these mechanisms can be combined in one term and written as Gt, where G is obviously the 
gradient of the net field and r is a coordinate point in the sample of a particular nucleus. In 
general, especially in a flow system, G and r will vary with time. In addition, there are statistical 
fluctuations throughout the system due to a myriad of phenomena. 

B = B 0 +B 1 +[G(r,t)T(t)]k (2.5) 
The precession frequency is proportional to the magnitude of the homogeneous magnetic 
induction B 0 . 

°>p=Yn|B 0 |=Y n B 0 (2.6) 

This simple relationship is the basis for a convenient technique, which is fruitful in obtaining 
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essential MRI data. The technique is to transform the observations from the laboratory reference 
frame into a rotating reference frame that rotates in such a way as to view the magnetization, M, 
as stationary. This is accomplished by constructing a rotating coordinate system (sometimes 
referred to as the body system and is usually denoted with primes on all coordinates to 
distinguish it from the lab frame) whose z'-axis is aligned with B 0 and whose x'-y' plane rotates 
with the precession frequency co p . It is as though B 0 has no effect on the magnetization. The 
advantage of observation from the perspective of the rotating system is that now all the necessary 
information needed to form an image is contained in the motion of the x', y' coordinates of the 
magnetization. That is the motion that will be caused by deviations from B 0 as a result of "extra" 
gradient fields — unwanted deviations, which produce image artifacts. 

In the body frame only the precession of M due to B] and the gradient term, Gt, will be 
observed. Bj will cause M to precess with frequency "Tn^i about the x'-axis in the y'z' 
plane. Thus in the rotating system the Bloch equation takes the form, 

m ^.(M.Bj-MiM.fc^i, (2 . 7) 

where k = k', and 

M^M^i' + My-j' + M^k', (2.8) 

Beff = B lx 4' + B 1/ j' + [G.r]k'. (2.9) 



Note that G*r is a scalar invariant. 

Since the intent is to work completely in the body system, all primes can be dropped — the 
rotating reference frame is implied. The Bloch equation and its component equations can now be 
written as 
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— - = y N (MxB eff ) ~ 

dt ±2 A i 



dM x ^ M x_ 
dt T 7 



dM. 



+ y N [G(r,t) • r(t)]M y - y N B ly M z , 



— = -Yn [G(r,t) • r(t)]M x - — - + Yn B 1x M z 
dt L 2 

dM , (M z -M 0 ) 

!2k = Yn (B ly M x - B lx M y ) • 

dt h 



(2.10) 



(2.H) 



(2.12) 



(2.13) 



This whole set of equations can be conveniently handled in matrix form. Defining two 
column vectors, |m) and |b), 



|m) = 



M x 
M y 
1VL 



(2.14) 



|B) = 



0 
0 

M 0 /T, 



(2.15) 



and the square matrix A 



A = 



- 1/T 2 Yn(Gt) -y N B ly 
-Y N (G.r) -1/T 2 Yn B 1x 
. Y N B ly -y n B 1x -1/T, 



(2.16) 



results in the compact matrix equation, 



d M > i \ i \ 

dt I / i / 



(2.17) 



This equation is the theoretical basis for most MRI measurements including measuring 
sequences of RF pulse measurements. 

MRI Imaging Equations With Flow 
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Whether it is patient movement, cardiac pulsing, or blood flow, motion of the target will 
cause inconsistencies in the phase and amplitude information obtained from measurements of the 
transverse magnetization, which can lead to image blurring and ghosts. It is possible to 
minimize the uncertainties in the data if the mechanisms that are producing them are known. 
However, in the case of blood flow, the flow properties themselves can be studied by visualizing 
the flow paths of arteries and veins (MR angiography) or by measuring the velocity of the fluid 
using direct flow measurements. There are essentially two approaches used for flow imaging. 
The first is called the phase-contrast method. It depends on the phase difference of the 
transverse magnetization of the flowing blood with respect to the stationary tissue surrounding 
the fluid. The second is called the modulus contrast method. This method is based on the 
difference of magnitudes of the transverse magnetization of the blood flow and the stationary 
tissue. 

Typical MRI systems involve echo-planar imaging: that is, a detection of spin echoes in 
planes or slices of the sample. An outline of the development of the basic echo-planar imaging 
equations, including flow, is given. Following their development, for illustrative purposes, the x 
and y components of Bi are set equal to zero, Bj x = Bi y = 0. The starting point for this 
development is with the definition of transverse magnetization given in terms of a complex 
coordinate system. This is done to facilitate a computational technique that makes use of Fourier 
transforms between real- and k-space. The transverse magnetization is thus defined as 

M T (t) = M x (t) + iM y (t) (3.1) 

where i = V-l • 

Meaningful data is obtained when the magnetic spins precess in phase. However, 
inhomogeneities in the magnetic field tend to de-phase the precessing spins. After the de- 
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phasing, an RF pulse, called a re-focussing pulse, is applied. The entire set of spins will begin to 
re-phase. At a time later, called the echo time, all spins precess with the same phase. The 
following analysis is based on echo spin detection. 

The total transverse magnetization of an excited slice, in which the gradient waveforms 
are functions of time, t, can be written as 



M T (t) = j|m(x,y)exp -ijVx,y;t')dt' 



dxdy, (3.2) 



where m(x, y) is the distribution of magnetization over the slice at a time just after excitation. 
And <j), which is a function, the "phase" function, which maps the time evolution of the 
transverse magnetization is given by 

♦(x,y;t) = y N [G(r,t)«r(t)] (3 j) 

Equation (3.2) assumes implicitly that there is no decay of magnetization due to spin-spin 
relaxation. For the case of no flow, r is essentially independent of time and the time integration 
over <[> becomes 

j4>(x,y;Odt' = k x x + k y y (3.4) 



where 



and 



then 



kx=Y N J 0 t Gx(Odt / , (3.5) 



k y=T N £G y (t')dt' (3.6) 



M T (t) = j|m(x,y;t)exp[-i(k x x + k y y)]dxdy = M T (k x ,k y ;t) (3.7) 



Equation (3.7) is the fundamental equation for MRJ. Since it has the form of a two-dimensional 
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Fourier transform an inverse Fourier transform can be defined: 

m(X ' y;t) = h JI M T( k x>k y ;t)exp[i(k x x + k y y)jdk x dk y 



(3.8) 



This transformation yields the required distribution of transverse magnetization for image 
construction in the selected slice at the time of the measurement. Spin-spin relaxation can now 
be easily incorporated by letting 

-ak x -Pk y 



M T (k x ,k y ;t)=>M T (k x ,k ;t)e-^ e 



(3.9) 



where a and |3 are decay constants inversely proportional to G X T 2 and G y T2, respectively. 

The incorporation of blood flow can easily be made by identifying a position point in the 
blood with r(f). Now r(t) can be expanded in a Taylor series about an arbitrary time coordinate, 
t e . The physical meaningful terms in the expansion are: 



r(t) = r(t e ) + (t-t e ) 



dr(t) 



dt 



+ «(t-t e )' 



,d 2 r(t) 



dt 2 



n3<* r(t) 
+ -( t " t c) -TT 2 



3!' 



(3.10) 



The first derivative corresponds to the velocity of the flow, the second derivative to the 
acceleration of the flow, and the third derivative to the rate of change of the acceleration of the 
flow, which would occur in turbulent flow. For illustrative purposes the acceleration is chosen 
as constant so that Eq. (3.10) reduces to the familiar kinematic equation of basic physics. 



r(t) = r(t e ) + v(t-t e ) + I a (t-t e ) 2 
Substituting Eq. (3.11) into (3.3) and using the definitions of Eqs. (3.5) and (3.6) yields 
J<|>(x,y;t')dt' » k • r e + k • v(t - t e ) + |k • a(t - t e ) 2 . 

For constant velocity, a = 0, the change in the integrated phase is just 

A<D= U(x,y;t')dt' - U(x,y;t')df =k«v(t-t e ). 

J Flow J NoFlow 



(3.11) 



(3.12) 



(3.13) 
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From this, flow parameters can be obtained from the information contained in the differential of 
the time integrated phase functions. The phase function is part of the general integral definition 
of the Fourier transform of the transverse magnetization of a slice given by Eq. (3.2). 

Bayesian Parameter Estimation Methods 
As a starting point, it is assumed that the MRI measurements result in a set of data that is 
represented by a column vector D of dimension 1 x N. If the data is multi-dimensional it is 
assume that the pixels, or voxels, can be ordered in such a way as to produce the data vector D, 
i.e., D(di,d2,.-*3 d>j). The data D can be modeled by a vector function, f, and n a noise 
component: 

D = f(a) + n, (4.1) 

a includes all the physical information, parameters, which arise in constructing a model to 
represent the data, n represents the system noise, which in general is assumed to be additive and 
Guassian in form. The model function many be determined from first principles if an adequate 
physical model of the system is already had, or phenomenologically from the system at hand. 
Indeed, it may also be determined by a complementary pairing of theoretical and empirical 
results. 

Because of the presence of noise, the inversion of Eq. (4.1) is non-unique and statistical 
procedures are necessary to obtain information about a. To this end the following probability 
densities are defined: 

P(a | DI) = Probability Density(confidence) that a is true, conditioned on D (4.2a) 
and any other prior information I. 

P(D | al) = Likelihood Function — Probability Density for the data D to have (4.2b) 
information concerning a. 
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P(a 1 1) = Prior Probability Density on a. (4.2c) 

P(D 1 1) - Probability Density for D. (4.2d) 
Note that, 

0<P(a|DI)<l, (4.3a) 

P(a | DI) =1 implies 100% confidence in a (4.3b) 

P(a | DI) = 0 implies a is ruled out (4.3c) 

An estimation of P(a |DI), also known as the "sampling distribution," is obtained by using Bayes 

theorem: 

P(a | DI) - [P(D | al) P(a 1 1)] / P(D 1 1) (4.4) 
For the case of Gaussian white noise the sampling distribution can be written as 

P(D | otl) = T] exp(-[D - f (a)] 1 2 -2 [D - f (a)]), (4.5) 



where r| is a normalization constant and S is the error covariant matrix for a. 

Two methods that have proven fairly successful and that are widely used in the 
estimation of a are the Maximum Likelihood (ML) method and the Maximum A Posteriori 
(MAP) method, which are special cases of the Bayesian statistical method. 

Maximum Likelihood (ML) Method 
This method corresponds to the case when no prior information is available for the 
estimation of a.. Then, the parameter vector a is estimated using 

V a P(D | al) = 0. (4.6) 

For Gaussian noise this becomes 

V a [D-f(a)] f 2- 2 [D-f(a)]-0. (4.7) 
That is Eq. (4.7) corresponds to the minimization of the nonlinear least-squares error. Thus one 
can use standard least-squares optimization procedures to estimate a. 
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Maximum A Posteriori (MAP) Method 
In this method it is assumed that a prior estimate for a, a 0 , is available and that the prior 
probability density is of Gaussian form. 

For the special case when both the likelihood and the prior estimate are Gaussian, a can 
be estimated from the following set of equations: 

P(a | DI) = a exp[-(a - a 0 ) T a a ~ 2 ( a - a 0 )], (4.8) 

and 

X - [D - f (a)] 1 T 2 [D - f (a)]+ (a - a 0 ) T a a " 2 ( a - a 0 ) (4.9) 
VaX = 0, (4.10) 
where a is a normalization constant and a a = Sy E, i.e., only the diagonal elements of Z. 

When the noise is taken as Gaussian the MAP method is equivalent to the minimization 
of %. Equation (4.10) results in an estimate of a clustered about a 0 . The clustering, or 
difference, will depend on just how close to the true value(s) a 0 was to begin with. (For the 
special case of Gaussian noise the MAP method turns out to be equivalent to Ridge Regession. 

Matched Filter Method 
The match filter method can generally be stated as the design of an "optimum filter + " for 
the received signal. (The filter is optimized based on the assumed best values of the 
predetermined filter parameters.) Here it is assumed that a sufficient, but not necessarily 
accurate, knowledge of the spin-spin relaxation time, T2, is available so that its affect may be 
incorporated into the transverse magnetization data via Eq. (3.9): 

MxCk^kyjO^MTCk^ky^e-^e^, (5.1) 

where, again, a and p are decay constants inversely proportional to G X T 2 and G y T 2 , respectively. 



18 



4683-107 

The factor of e~ ak * e~^ ky is called the "filter function" in this method. The values a and (3 are 
predetermined and fixed in the calculation, A two-dimensional power spectrum, P(k x ,k y ;t), may 
easily be obtained by taking the square of the absolute magnitude of the (complex) transverse 
magnetization, 



P(k x5 k y ;t): 



M T (k x ,k y ;t)e- ak »e- pk 12 



(5.2) 



The position of the peaks in the spectrum as a function of time can be used to obtain information 
about the velocity of the flow. 

The matched filter method requires a priori knowledge of the decay constants a and (5 
and the analysis critically depends on the choice of these parameters. In the Bayesian method, 
previously described, the optimal values of these parameters are determined as part of the 
calculation. The flexibility of being able to adjust the parameters in accord with the actual test 
data is a decided advantage over the use of predetermined, fixed parameters employed by current 
methods in use such as with the matched filter method. 

Under certain conditions it is advantageous to use a radially averaged power spectrum, 
defined by 

( P Me=^f P ( k ^t)<ie> (5.3) 

where 0 arises through the direction cosines defining k x = k Cos(9) and k y = k Sin(G) in the plane 
of the slice on which measurements are being made. Information concerning flow velocity is 
readily extracted from this procedure as well Compared to Eq. (5.2), the radially averaged 
spectrum is simpler to deal with since the number of variables is reduced. 

Bayesian Method 
The fundamental equation of MRI given in Eq. (3.7) 
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M T (k x ,k y ;t) = JJm(x,y;t)exp[-i(k x x + k y y)]dxdy . (5.4) 

For illustrative purposes it is assumed that the object function, the inverse Fourier transform of 
Eq. (5.4), takes the form of a multiple Dirac delta function distribution, 

m(x,y;t)~ A8(x - x 0 )8(y - y 0 ), (5.5) 

where 

Aocm(x 0 ,y 0 ;t). (5.6) 
For this special case Eq. (5.4) reduces to 

___Mt^ Aexp[-i(k x x 0 + k y y 0 )]exp(-ak x )exp(-pk y ), (5.7) 
where A, x 0 , y 0 , a and p are the parameters to be estimated. MRI data, D 5 is then modeled as 

D(k x ,k y ;t) = M T (k x , k y ;t) + n(k x ,k y ), ( 5 . 8 ) 
where r\ is the noise function and M T is given by Eq. (5.1). An estimation of these parameters 
follows the general procedures, via Bayesian Probability Theory. The techniques for static 
imaging can now be applied to Eq. (5.8) for MRI flow imaging. Also, estimation of a and p via 
the Bayesian approach will return highly probable values and thus an accurate estimate of T 2 — 
which is a characteristic of the patient's system under study. x 0 and y 0 contain flow velocity 
information and an accurate estimation of them, via the Bayesian approach, is critical for an 
evaluation of the patient's condition. 

The general formulation of Magnetic Resonance Imaging was constructed with the use of 
the Bloch equations. The Bloch equations are critical for the analysis because they account for 
the dephasing of the nuclear magnetic spin system due to spin-spin and spin-lattice interactions. 
In-phase spin precessional motion is essential to obtain meaningful data— an estimate of the 
amount of dephasing is a necessity for an accurate interpretation of the data. The effects of spin- 
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lattice and spin-spin interactions were accounted for in the Bloch equations by their associated 
relaxation, or decay, times, Ti and T2, respectively. 

Prior art development of MRI was for stationary targets. That development has been 
extended to incorporate the motion of flowing targets such as blood in a vessel. Typically, the 
imaging is done in a process called echo-planar imaging in which data is simultaneously taken in 
multiple slices of the target system. A convenient technique for analyzing the data is to work in 
reciprocal or Fourier space. Once the transformation was made, the incorporation of flow was 
made by multiplying the function by exponential functions decaying in Fourier space at a rate 
proportional to the spin-spin relaxation rate — Eq. (3.9). To complete the picture of flow a 
kinematic equation describing the motion of a position vector within each slice was derived. 
This was then related to the "phase" function given as part of the general integral definition of 
the Fourier transform of the transverse magnetization of a slice given by Eq. (3.2). It was then 
shown that the differential of the time integrated phase functions, with and with flow included, 
contains all the necessary flow parameter information. 

The Bayesian methodology used for parameter estimation of any problem in which 
critical parameters are to be extracted from a set of given data was presented. The Bayesian 
method is unique among all methods, which strive to estimate such parameters in that it is a 
dynamical method based on conditional probabilities. As a result the parameters can be adjusted 
for the best fit of the data in real time with minimum uncertainty. Moreover, their values may be 
improved upon by the addition of a new knowledge of the physical situation in the form of 
empirical and/or theoretical models that may include features such as boundary conditions or 
perturbations. This is all accomplished by inputting a model function and the use of conditional 
probabilities based on Bayes' Theorem. As an illustration, the Bayesian approach was applied to 
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the Method of Maximum Likelihood and the Maximum A Posteriori (MAP) Method. 

Lastly, the material was unified and applied to the problem MRI flow parameter 
estimates. For comparison, an example of the standard technique using matched filters was 
outlined. Clearly the starting place for both the matched filter method and the Bayesian method 
was with the same model function, Eq. (5,1). However, the matched filter method requires a 
priori knowledge of the decay constants a and p and the analysis critically depends on the 
choice of these parameters. Essentially, the analysis can, and often does, get carried out with 
only a sufficient, but not necessarily accurate, knowledge of the input parameters — in this case 
those related to the spin-spin relaxation time. 

In the Bayesian method the optimal values of the input parameters as well as the 
parameter of interest are determined as part of the calculation. The input values are used only as 
an initial guess to start the algorithm. The flexibility of being able to adjust the parameters in 
accord with the actual experimental data is a decided advantage over the use of predetermined, 
fixed parameters employed by current methods in use such as with the matched filter method. 
Moreover, as was pointed out earlier the values of the parameters may be improved upon by the 
addition of a new knowledge of the physical situation in the form of empirical and/or theoretical 
models that may include features such as boundary conditions or perturbations. The Bayesian 
method is a dynamic and unique method for parameter estimation, and prior to this disclosure 
has not been applied to MRI flow measurements. A summary flow chart of the methodology just 
described is shown in Fig.2. MRI data 200, an MRI model function f(a) and prior information 
204 are processed utilizing Bayesian methodology 206 to image the feature of interest a. This 
feature of interest a is utilized in a clinical application 210 to make a medical decision 212. 



22 



4683-107 

Unlike traditional methods, which assume an underlying Gaussian noise for the data the 
Bayesian methods do not assume a specific noise model. Rather the Bayesian method examines 
the data to infer the relevant noise model appropriate for the problem. It also will compare the 
probabilities for several noise models and estimate, which noise model one is validated by the 
data. Thus the Bayesian method can outperform traditional methods in this respect if the 
underlying noise model is non-Gaussian and nonstationary. In a sense the data speaks for itself 
through the Bayesian method. 

In view of the foregoing description, numerous modifications and alternative 
embodiments of the invention will be apparent to those skilled in the art. Accordingly, this 
description is to be construed as illustrative only and is for the purpose of teaching those skilled 
in the art the best mode of carrying out the invention. Details of the structure may be varied 
substantially without departing from the spirit of the invention, and the exclusive use of all 
modifications, which come within the scope of the appended claim is reserved. 
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